function Sigma= NW_hac(vars)
%calculates Newey-West covariance matrix based on number lags set by "lags"
%below
lags=4;

Sigma0=cov(vars);
Sigma_cov=@(k) (1/(size(vars,1)-k))*(vars(1:end-k,:)-repmat(mean(vars(1:end-k,:)),size(vars(1:end-k,:),1),1))'*(vars(1+k:end,:)-repmat(mean(vars(1+k:end,:)),size(vars(1+k:end,:),1),1));

Sigma=Sigma0;
for n=1:lags
    Sigma=Sigma+(1-n/(lags+1))*(Sigma_cov(n)+Sigma_cov(n)');
end
end

